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Presently,  Rotational  Isomeric  State  (RIS)  theory  directly  addresses  polymer  chain  conformation  as 
it  relates  to  mechanical  response  trends.  The  primary  goal  of  this  work  is  to  explore  the  adaptation 
of  this  methodology  to  the  prediction  of  material  stiffness.  This  multi-scale  modeling  approach 
relies  on  ionomer  chain  conformation  and  polymer  morphology  and  thus  has  potential  as  both 
a  predictive  modeling  tool  and  a  synthesis  guide.  The  Mark-Curro  Monte  Carlo  methodology  is 
applied  to  generate  a  statistically  valid  number  of  end-to-end  chain  lengths  via  RIS  theory  for  four 
solvated  Nation  cases.  For  each  case,  a  probability  density  function  for  chain  length  is  estimated 
using  various  statistical  techniques,  including  the  classically  applied  cubic  spline  approach.  It  is 
found  that  the  stiffness  prediction  is  sensitive  to  the  fitting  strategy.  The  significance  of  various 
fitting  strategies,  as  they  relate  to  the  physical  structure  of  the  polymer,  are  explored  so  that  a 
method  suitable  for  stiffness  prediction  may  be  identified. 


1.  Introduction 

Ionic  polymers  comprise  the  active  layer  in  Ionic  Polymer-Metal  Composites  (IPMCs),  which  were 
first  identified  just  over  a  decade  ago  [1-5]  and  now  constitute  an  emerging  class  of  soft  transducers 
(see  Figure  1).  Because  they  generate  large  strain  in  response  to  low  electric  field  stimulation  and 
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Figure  1:  Cross-sectional  view  of  a  typical  IPMC  with  illustration  of  the  assumed  ionic  polymer 
clustering  morphology  based  on  Hsu  and  Geirke  [12]. 

have  high  gravimetric  energy  density,  there  has  been  considerable  conjecture  over  the  potential 
applications  of  these  soft  transducers  [2,  3,  6-8]  and  considerable  fruitful  investigations  towards 
these  ends  [9-11].  The  high  gravimetric  energy  density  is  a  function  of  both  the  electromechanical 
transduction  mechanisms  of  the  ionic  polymer  layer  and  the  global  IPMC  stiffness.  Most  models  of 
the  transduction  behavior  rely  on  the  cluster  morphology  of  the  ionic  polymer  layer,  as  proposed  by 
Hsu  and  Gierke  [12],  where  the  material  has  been  solvated  with  water.  In  brief,  the  backbone  of  the 
ionic  polymer  chain  is  hydrophobic  whereas  the  side  chains  terminate  in  hydrophilic  ionic  groups. 
Hsu  and  Geirke  propose  a  clustering  of  these  hydrophilic  ionic  side  groups  and  the  water  that 
has  been  taken  up  by  the  material.  The  model  further  suggests  an  idealized  structuring  whereby 
the  clusters  are  of  essentially  constant  radius,  uniformly  distributed  throughout  the  material,  and 
interconnected  by  channels.  Subsequently,  many  have  sought  to  quantify  the  cluster  radius  and 
spacing  as  functions  of  polymer  equivalent  weight,  hydration  level,  counterion,  etc.  —  see  for 
instance  the  modeling  works  of  Datye  and  coworkers  [13, 14]  and  the  experimental  work  of  Lehmani 
et  al.  [15]  and  Barbi  et  al.  [16]. 

Current  transduction  models  tend  to  be  categorized  in  two  classes:  models  assuming  a  hydraulic 
description  of  the  electromechanical  response  [17,18]  and  hybrid  models  including  electrostatic,  elas¬ 
tic,  and  osmotic  effects  [8, 19,20].  Despite  the  disagreement  over  the  underlying  electromechanical 
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mechanisms,  the  resulting  predictions  of  all  of  the  above  modeling  approaches  are  dependent  upon 
ionomer  stiffness.  Moreover,  it  has  been  noted  that  the  ionomer  stiffness  affects  transducer  re¬ 
sponse  both  directly,  via  classic  energy  density  consideration,  as  well  as  indirectly,  via  an  impact 
on  transport  properties  [20,21].  Thus  a  modeling  method  which  can  predict  ionomer  stiffness  as  a 
function  of  its  specific  composition  (ionomer  type,  counter  ion  type,  solvation  type  and  level,  etc.) 
is  desirable. 

Existing  stiffness  prediction  approaches  fall  short  of  this  objective  by  being  either  too  case 
specific  or  too  general.  For  instance,  as  described  by  Treloar  [22],  the  variation  in  rubber  stiffness 
as  a  function  of  fluid  uptake  has  long  been  estimated  by 

Gwet  =  GdryVl/\  (1) 


where  G  is  the  shear  modulus  and  v  is  the  volume  fraction  of  rubber  in  the  rubber-liquid  mix. 
Further,  incompressibility  is  a  classic  assumption  of  rubber  elasticity;  thus  Poisson’s  ratio  is  set  to 
0.5  and  Young’s  modulus  may  also  be  directly  solved.  Unfortunately,  experimental  data  indicates 
that  while  many  aspects  of  ionomer  mechanical  behavior  parallel  that  of  rubber,  the  way  in  which 
solvent  is  taken  up  (by  way  of  clustering)  is  fundamentally  different  than  for  a  classic  rubber,  and 
hence  the  above  general  relation  becomes  invalid.  A  second,  empirical  approach  is  proposed  by  Li 
and  Nemat-Nasser  [19],  based  on  the  works  of  Hsu  and  Geirke  [12]  and  Grot  et  al.  [23].  For  the 
specific  case  of  Nafion,  Young’s  modulus  is  given  by 


E(cw)  =  E0  exp 


/  100cwpw  1200  -EWV 
V(l-cw)pd  +  20  )\  ’ 


(2) 


where  Eq  =  275  MPa,  a  =  0.0294,  cw  is  the  water  volume  fraction,  pw  and  pd  are  the  densities 
of  water  and  dry  polymer,  and  EW  is  the  equivalent  weight  of  the  Nafion  sample  in  question. 
Example  predictions  for  this  empirical  relation  are  131  MPa  and  100  MPa  for  fully  hydrated  (via 
boiling)  Nafion  1100  equivalent  weight  with  sodium  and  lithium  as  the  counterions,  respectively. 
When  compared  to  the  experimentally  measured  values  of  80  MPa  and  75  MPa  respectively  [8], 
this  estimate  is  good  for  the  current  state  of  the  art;  however,  it  has  been  formulated  specifically 
for  Nafion  rather  than  for  ionomers  in  general. 
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Rotational  Isomeric  State  (RIS)  theory  as  described  by  Flory  [24]  has  long  been  used  as  a  means 
to  anticipate  mechanical  response  trends  in  polymeric  materials  based  on  the  underlying  polymer 
chain  conformation.  The  fundamental  idea  of  this  approach  is  that  any  given  bond  within  a  single 
polymer  chain  is  limited  to  a  discrete  number  of  possible  bond  angles,  based  on  corresponding  low 
energy  states.  From  this  base,  any  number  of  statistical  methodologies  may  be  applied  to  anticipate 
final  chain  conformation  and  corresponding  response  to  load.  As  discussed  by  Curro  and  Mark  [25], 
application  of  the  Gaussian  approach  to  RIS  theory,  as  described  by  Flory,  does  not  represent  a 
reasonable  approximation  for  short  chain  polymers  displaying  rubber  elasticity.  A  more  realistic 
approach  is  to  apply  RIS  theory  in  combination  with  a  Monte  Carlo  methodology,  from  which  a 
large  number  of  end-to-end  distances  r  are  generated.  The  r  values  can  then  be  used  to  estimate  an 
appropriate  probability  density  function  P(r)  for  the  short  chain  case.  A  cubic  spline  approximate 
provides  one  technique  that  has  been  employed  to  estimate  P(r). 

The  Mark-Curro  approach  has  also  been  applied  to  investigations  on  the  effect  of  particle 
reinforcement  [26-28].  In  these  particle  reinforcement  investigations,  various  volume  fractions, 
inclusion  sizes,  and  inclusion  shapes  and  orientations  are  considered,  where  the  chain  conformation 
is  excluded  from  occupying  any  volume  dedicated  to  an  inclusion.  Thus,  the  predicted  constitutive 
responses  are  due  solely  to  the  altered  conformations  of  the  polymer  matrix  and  do  not  account  for 
the  load  carrying  capacity  of  the  inclusions.  A  later  preprint  [29]  suggests  that  efforts  are  underway 
to  account  for  the  load  carrying  capacity  of  the  inclusions  as  well,  and  subsequently  the  methodology 
may  be  used  for  stiffness  predictions.  To  date  this  method  has  not  been  used  for  stiffness  predictions 
for  reasons  largely  related  to  the  limitations  introduced  by  way  of  the  simplifying  assumptions  (i.e., 
the  commonly  applied  3-state  models  of  chain  bond  angles  are  unable  to  capture  helical  chain 
coiling  [30];  uncompensated  free-end  effects,  accurate  identification  of  polymer  molecular  weight, 
etc.  [22]). 

In  this  paper  it  is  hypothesized  that  reasonable  estimations  for  the  above  noted  limitations  are 
possible  so  long  as  a  more  appropriate  statistical  approach  is  first  employed.  Stiffness  predictions 
are  thus  considered  for  various  similar,  solvated  ionomers. 

It  is  understood  that  ionomer  stiffness  is  a  function  of  multiple  parameters  including  ionomer 
type,  solvation  type  and  level,  and  counterion  type.  In  this  work  RIS  theory  is  applied  in  a  manner 
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analogous  to  the  particle  inclusion  works  [26-28]  in  order  to  capture  each  of  these  contributions, 
where  inclusion  volumes  are  now  taken  to  be  clusters.  Cluster  morphology  is  assumed  a  priori, 
based  on  available  experimental  measurements  for  the  various  cases;  the  subsequent  backbone 
conformation  in  response  to  this  morphology  is  applied  to  obtain  stiffness  predictions.  In  addition, 
it  is  assumed  that  these  clusters,  under  specified  conditions,  may  act  as  backbone  cross-linking 
junctions,  and  thus  the  ionomer  stiffness  is  based  on  polymer  short  chain  response.  Moreover, 
this  computational  approach  permits  direct  elimination  of  free-end  effects  and  provides  a  means  to 
estimate  a  representative  cluster-to-cluster  molecular  weight  (or  effective  molecular  weight  based 
on  the  distance  between  cluster  communication  points  on  the  polymer  backbone). 

The  significant  differences  between  the  particle  reinforcement  investigations  and  the  current 
technique  are  the  following:  (i)  the  inclusions  are  hydrophilic  clusters  with  which  the  hydrophobic 
backbone  may  interact  via  the  pendant  chains,  thus  affecting  the  end-to-end  distances  considered, 
(ii)  the  probability  density  function  P(r)  is  estimated  using  various  statistical  approaches,  and  (iii) 
the  various  predictions  are  considered  both  in  terms  of  the  predicted  trends  and  the  magnitude  of 
the  predicted  stiffness.  It  is  found  that  stiffness  predictions  are  sensitive  to  the  choice  of  statistical 
methodology  and  thus  represent  a  tool  for  expanding  the  capabilities  of  this  modeling  approach 
to  stiffness  predictions.  Details  of  the  various  statistical  approaches  are  discussed  in  [31].  In  this 
work,  the  focus  is  on  the  impact  these  various  methodologies  have  on  the  predicted  response. 

The  paper  is  organized  as  follows.  The  second  section  discusses  the  RIS  model  construction  for 
the  ionic  polymer  Nafion  in  four  hydrated  forms:  1100  Equivalent  Weight,  in  Na+  and  Li+  forms, 
with  two  assumed  cluster  distributions  each.  Section  three  briefly  describes  the  statistical  analysis 
approaches  and  the  fourth  section  presents  model  predictions  and  discussion. 

2.  Multiscale  Model  Development 

The  mechanical  response  of  two  hydrated  Nafion  1100  Equivalent  Weight  cases  in  salt  form  are 
considered,  one  with  Na+  and  the  other  with  Li+  counterions,  where  it  is  assumed  that  the  ion  ex¬ 
change  is  complete.  The  Mark-Curro  approach  to  RIS  theory  [25]  is  applied  in  a  manner  analogous 
to  that  employed  in  investigations  focused  on  on  the  effect  of  particle  reinforcement  [26-28].  In  the 
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Figure  2:  Chemical  composition  of  a  Nation  monomer. 

particle  reinforcement  studies  the  chain  conformation  is  excluded  from  occupying  any  volume  ded¬ 
icated  to  a  rigid  inclusion.  In  this  investigation  the  inclusions  take  the  form  of  hydrophilic  clusters 
rather  than  rigid  particles.  The  significant  difference  in  the  nature  of  the  inclusion  is  addressed 
by  way  of  the  additional  constraint  that  the  backbone  may  have  discrete  points  of  communication 
with  the  hydrophilic  clusters  via  the  pendant  chain  attachment  points.  These  attachment  points 
are  assumed  to  approximate  cross-linking  points,  and  thus  mark  an  end-point  for  the  given  chain 
length  r. 

2.1.  Simulation  Model  for  Estimating  Chain  Length 

To  implement  this  model,  appropriate  statistical  weight  matrices  for  backbone  chain  conforma¬ 
tion,  bond  lengths,  and  distances  between  pendant  chain  attachment  points  must  be  identified. 
The  repeating  monomer  unit  in  the  chemical  structure  of  Nation  is  depicted  in  Figure  2,  where 
(CF2CF)(CF2CF2)n  represents  the  repeating  unit  of  the  backbone  and  SOJ  is  the  pendant  chain 
terminal  ionic  group  which  tends  toward  the  hydrophilic  clusters.  If  the  87/13  Nafion  case  is  con¬ 
sidered,  meaning  that  there  are  approximately  13  (CF2CF)  groups  to  every  87  (CF2CF2)  groups 
in  the  total  length  of  the  backbone,  then  n  is  most  often  approximately  7.  In  other  words,  the 
bulk  of  the  backbone  takes  the  form  of  PTFE,  or  Teflon  [32,33].  Because  it  is  the  backbone  of  the 
polymer  which  is  responsible  for  the  load  carrying  capacity  of  the  ionomer,  the  PTFE  conformation 
is  assumed  for  all  bond  placements  except  pendant  chain  attachment  points. 

Whereas  a  typical  value  for  n  in  the  chemical  structure  of  Nafion  is  approximately  7,  there  is  in 
fact  variability  in  this  value.  To  account  for  this  variation,  n  is  sampled  from  a  discrete  probability 
distribution  with  a  minimum  of  5,  a  maximum  of  11,  and  a  mean  of  approximately  7  [31]. 

Still  to  be  addressed  is  the  number  of  times  the  repeat  unit  of  Figure  2  is  expected  for  a 
single  polymer  strand  of  Nafion.  Because  the  actual  molecular  weight  is  not  readily  available, 
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Figure  3:  Depiction  of  the  in-plane  bond  angle  6  =  116°  and  PTFE  backbone  component  length 
lc~c  =  1-53  A.  The  angle  6'  is  the  supplementary  angle  to  6  having  the  value  64°. 

in  part  because  polymers  do  not  have  fixed  chain  lengths,  a  distribution  of  lengths  is  assumed. 
The  assumed  distribution  anticipates  that  total  chain  lengths  contain  between  135  and  225  repeat 
units.  Note,  however,  that  the  lack  of  available  information  on  this  parameter  is  considered  non- 
critical  because  this  material’s  load  bearing  capability  is  based  on  the  distance  between  cluster 
communication  points  along  a  given  backbone  rather  than  on  the  entire  backbone  chain  length. 

The  model  progression  begins  with  the  placement  of  a  cluster  configuration  in  a  3-dimensional 
grid  of  (5000  A)3.  The  grid  size  is  selected  to  accommodate  the  longest  anticipated,  fully  extended 
backbone.  Two  cluster  distribution  schemes  are  considered  as  detailed  in  the  preliminary  work  [31]. 
The  first  relies  directly  on  the  simplified  model  of  Hsu  and  Gierke  [12]  where  the  clusters  are 
uniformly  distributed  on  a  square  grid  analogous  to  the  cubic  crystalline  structure.  The  second 
uniformly  distributes  clusters  of  the  same  radius  and  net  volume  fraction  as  the  first  case,  but  on 
a  staggered  grid  analogous  to  the  HCP  crystalline  structure.  The  cluster  size  and  volume  fraction 
for  the  cases  studied  are  provided  in  Table  1.  For  simplicity,  calculation  of  the  volume  taken  up  by 
the  ionic  groups  is  currently  neglected.  Should  this  prove  significant,  the  estimation  proposed  by 
Hsu  and  Gierke  [12]  can  be  applied  to  cluster  redistribution. 

Once  the  cluster  distribution  is  assigned,  a  single  PTFE  bond  is  randomly  placed  in  the  3- 


Table  1:  Cluster  Morphology  Values  based  on  [8]. 


Li+ 

Na+ 

Cluster  Volume  Fraction 

38% 

30% 

Cluster  Radius 

23  A 

21  A 

Cluster  Center-to-Center  Separation 

50.6  A 

50.4  A 
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dimensional  grid,  where  cluster  locations  are  excluded.  For  PTFE  the  in-plane  bond  angle  and 


length  are  fixed  at  6  =  116°,  and  lc~c  =  1-53  A,  respectively,  as  illustrated  in  Figure  3.  Thus, 
the  in-plane  angle  of  the  second  bond  is  known;  however,  a  statistical  weight  matrix  must  be 
applied  in  order  to  identify  an  appropriate  out-of-plane  rotation  angle  ,  where  it  is  understood  that 
various  angular  orientations  between  bonds  correspond  with  local  low  energy  conformations.  Each 
subsequent  bond  must  be  similarly  placed.  For  PTFE,  both  3-state  and  4-state  statistical  weight 
matrices  have  previously  been  studied,  where  the  4-state  matrices  have  been  shown  to  be  somewhat 
preferred  [30].  For  the  4-state  case  the  possible  out-of-plane  angles  are  +15°, +120°, —120°,  and 
—15°  (or  trans+,  gauche+,  gauche-, trans-).  The  4-state  matrices  applied  are 
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where  a  =  0.2,  a'  =  2.0,  c v  =  0.2,  f3  =  0.5  at  25°C  [30],  and  the  rows  and  columns  are  indexed 
in  the  order  +15°,  +120°,  —120°,  and  —15°.  After  these  matrices  have  been  normalized  [31],  the 
orientation  of  the  second  bond  is  dictated  by  C/2,  the  third  by  C/3,  the  last  by  C/jv,  and  all  others 
by  C/fc,  with  two  exceptions.  The  matrices  are  overridden  when  a  bond  placement  coincides  with  a 
cluster  or  when  the  bond  corresponding  to  a  pendant  chain  attachment  point  is  being  placed.  The 
first  exception  addresses  the  reality  that  the  backbone  chain  is  hydrophobic  and  will  reconfigure 
upon  solvation  to  an  alternate  low  energy  state  in  order  to  avoid  the  hydrophilic  region.  The  second 
exception  addresses  the  reality  that  the  terminal  ionic  group  will  tend  toward  a  hydrophilic  region 
and  again  locally  reconfigure  the  chain  upon  solvation.  It  is  understood  that  application  of  the 
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Figure  4:  Distance  r  between  cluster  interaction  points  for  a  single  Nation  polymer  chain. 

in-plane  PTFE  bond  angle  in  this  case  is  a  simplifying  assumption. 

When  a  bond  placement  would  locate  the  next  bond  within  a  cluster,  another  random  number  is 
generated  to  find  an  alternate  but  viable  low  energy  conformation.  If  all  out-of-plane  rotation  angles 
result  in  cluster  coincidence,  the  chain  is  terminated.  When  a  bond  corresponding  to  a  pendant 
chain  connection  point  is  being  placed,  of  the  possible  ±15°  and  ±120°  out-of-plane  PTFE  rotations, 
the  angle  that  minimizes  the  distance  between  the  carbon  atom  with  the  attached  side  chain  and 
the  nearest  cluster  is  used.  It  is  assumed  that  the  pendant  chain  terminal  ionic  group  successfully 
communicates  with  the  nearest  cluster  if  the  distance  between  the  backbone  connection  point  and 
the  nearest  cluster  is  within  8  A.  This  value  is  an  estimate  of  the  fully  extended  length  of  the 
pendant  chain.  It  is  found  that  application  of  this  value  leads  to  between  75%  and  92%  successful 
cluster  communication  depending  on  the  assumed  cluster  configuration.  It  is  not  expected  that 
100%  of  the  ionic  groups  will  reside  within  a  cluster  and  thus,  this  range  is  deemed  reasonable. 

Each  simulated  Nafion  backbone  chain  can  have  multiple  r  values  as  illustrated  in  Figure  4. 
Furthermore,  the  first  and  last  r  value  for  each  chain,  or  free  ends,  do  not  support  any  load,  and 
therefore  do  not  contribute  to  stiffness.  Predictions  with  and  without  correction  of  the  free-end 
effects  are  considered. 

In  order  to  assure  statistical  validity,  a  large  number  of  r  values  (about  10,000)  are  generated. 
From  this,  an  expression  for  the  chain  length  probability  density  function  P(r)  must  be  estimated. 
Development  of  an  appropriate  expression  is  the  subject  of  Section  3. 

2.2.  Macroscopic  Model  for  Estimating  Stiffness 

Assuming  an  appropriate  P(r)  expression  has  been  developed,  per  Boltzmann’s  approach  to  sta¬ 
tistical  thermodynamics,  it  can  be  directly  related  to  the  entropy  of  the  chain  as  a  function  of  the 
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number  of  conformations  available  for  a  given  state  [25] 


S(r)  =  c  +  k\nP(r) 


(4) 


where  k  is  Boltzmann’s  constant  and  c  is  a  constant  of  integration  which  drops  out  when  the 
difference  in  entropy  for  the  unperturbed  configuration  is  taken  with  respect  to  the  distorted  con¬ 
figuration.  In  this  approach  it  is  assumed  that  under  load  the  rotation  about  the  individual  bonds 
is  unrestricted  which  allows  the  assumption  that  the  Helmholz  free  energy  is  strictly  a  function 
of  entropy.  For  the  assumptions  of  rubberlike  elasticity,  the  “three  chain”  model,  as  described  by 
Treloar  [32],  yields  the  relation 


AS 


S(rQa)  +  2 S(rQa  1/2)  -  3 S(rQ)  , 


(5) 


for  the  change  in  entropy  upon  distortion.  Here  v  is  the  number  density  of  network  chains,  rQ  is  the 
root  mean  square  of  the  r  values  applied  in  the  development  of  P(r),  and  a  =  L/ Li  is  the  relative 
length  of  the  sample.  The  nominal  stress  f*  is  then  given  by 


r  =  -T 


f  dAS\  _  vkTr0 
\  da  )  T  3 


G'{r0a)  —  a  3^2G'(r0a  1/2) 


where  G{r)  =  lnP(r)  and  G'(r)  =  - .  The  corresponding  modulus  [/*]  is  given  as 


[/*] 


f* 

a  —  a~2 


(6) 


(7) 


For  small  strains  (a  — >  1),  the  modulus  in  (7)  approaches  Young’s  modulus  for  the  polymer  matrix 
alone;  this  small  strain  modulus  will  be  referred  to  as  iqnjt.  The  global  Youngs’s  modulus,  Pave, 
is  expected  to  be  lower  than  E-m it  due  to  the  volume  fraction  of  hydrophilic  clusters  which  are  not 
expected  to  be  load  bearing.  Because  only  the  global  stiffness,  F'avej  is  experimentally  available  the 
means  to  properly  express  -Emit,  and  subsequently  relate  it  to  Eave  are  required.  Consider  first  the 
expression  for  E\n\t. 
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Substituting  equation  (6)  into  equation  (7)  and  applying  L’HopitaPs  rule,  yields  the  relation 


r  _  ukTr0  \  r0P{r0)P"(r0)  -  rQ  [ P'{r0 )]2  +  P{r0)P'{r0) 

6  1  [Pir o)]2 

for  the  initial  stiffness. 

In  order  to  relate  to  E&ve,  it  is  first  noted  that  E-m\t  is  directly  analogous  to  matrix  stiffness, 
Em,  of  classic  composite  materials  theory.  Thus,  a  modified  Voigt-Reuss  approximation  for  plane 
stress  [34]  is  used  to  relate  £init  to  the  experimentally  available  values,  E ave.  This  approximation 
is  given  by 


E \ 


init 


—  E  — 

—  — 


-^ave  — 


(9) 


where  /  is  the  volume  fraction  not  taken  up  by  water,  which  can  be  equivalently  expressed  in  terms 
of  the  normalized  volume  after  uptake  Vp  (i.e.,  if  upon  hydration  the  volume  increases  by  44%,  then 
Vp  is  1.44). 

Traditionally  RIS  has  been  used  for  the  purpose  of  studying  mechanical  response  trends  where 
model  predictions  described  by  equations  (6)  and  (7)  are  normalized  with  respect  to  vkT.  However, 
a  value  must  be  assigned  to  this  parameter  in  order  to  predict  stiffness.  Identification  of  the  number 
density  of  chains  v  requires  the  application  of  some  simplifying  assumptions.  To  begin,  v  is  given 
by 


pN A 

"  “  MWC_C 

where  p  is  the  material  density,  N\  is  Avagadro’s  number,  and  MWC_C  is  the  molecular  weight  of 
the  portion  of  the  chain  between  two  cluster  communication  points. 

The  material  density  is  a  readily  available,  experimentally  determined  value,  whereas  Avagadro’s 
number  is  a  constant.  This  leaves  MWC_C,  which  is  estimated  based  on  the  percentage  of  successful 
pendant  chain-to-cluster  communications.  The  estimate  assumes  that  n  =  7  in  the  bulk  of  the 
repeat  units  (Figure  2),  which  corresponds  to  a  molecular  weight  of  that  repeat  unit  of  781  g/rnol, 
where  again,  only  the  backbone  is  considered  to  be  load  bearing.  It  is  next  assumed  that  failed 
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cluster  communications  are  sufficiently  infrequent  that  it  is  highly  improbable  that  there  would 
be  two  failures  in  sequence.  Thus,  a  single  failure  corresponds  to  a  local  MWC-C  of  1562  g/rnol. 
Applying  these  values  as  a  ratio  of  their  occurrence  is  then  used  to  estimate  MWC_C.  For  example, 
if  the  failure  rate  is  15%,  then  MWC_C  ~  0.85(781)  +  0.15(1562)  =  898  g/mol. 


3.  Statistical  Analysis 


The  details  of  the  statistical  fitting  methods  applied  to  the  development  of  various  P(r)  expressions 
for  the  same  distribution  of  r  values  are  discussed  in  [31].  Provided  here  is  a  brief  overview  of  the 
cubic  spline  and  Johnson  family  methods.  In  the  traditional  approach  of  RIS  Monte  Carlo  theory, 
a  cubic  spline  approach  is  used  to  estimate  P(r).  The  simulation-generated  r  values  are  divided 
into  groups,  or  bins,  and  a  knot  value  is  specified  for  each  bin  so  that  on  each  interval  (delimited 
by  the  knots),  a  cubic  polynomial  is  fit  to  the  data.  If  we  define  the  indicator  function  L*(?’ )  for 
the  subinterval  (fcj,fc*+ 1],  where  ki  is  the  itli  knot,  as 


U{r) 


1,  for  ki  <  r  <  ki+i, 

< 

0,  otherwise, 


(11) 


then  the  result  is  a  piecewise-polynomial  estimate  of  P(r)  of  the  form 

1  n* 

P (r)  =  ^  E  L *(r)  iA^r  -  k^  +  B^r  ~  ki )2  +  A(r  ~  k i)  +  A]  ,  (12) 

i= 1 

where  n*  is  the  number  of  knots,  A,  B,  C,  and  D  are  the  corresponding  coefficients  of  the  cubic 
equations,  and  K  is  a  normalizing  constant  so  that  the  function  P(r)  given  in  (12)  integrates 
to  unity.  While  this  expression  is  easily  differentiable  for  substitution  into  equations  (4)-(7),  one 
challenge  with  this  approach  is  that  cubic  polynomials  are  often  not  sufficiently  flexible  to  accurately 
model  the  tail  behavior  of  the  density.  Additionally,  the  shape  of  the  estimated  density  is  highly 
dependent  on  the  number  and  location  of  the  knots.  For  example,  if  the  number  of  knots  is  large, 
the  corresponding  estimate  of  P(r)  not  only  gives  the  visual  impression  of  being  unnaturally  noisy, 
but  also  yields  nonphysical  response  trends.  Finally,  the  details  of  the  fitting  methodology  (such 
as  knot  selection)  are  typically  not  well-defined.  In  the  current  work,  a  series  of  fits  are  considered 
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where  the  number  and  location  of  knots  is  varied,  in  order  to  illustrate  the  physical  implications  of 
using  a  cubic  spline  approach. 

Alternatively,  the  application  of  the  Johnson  family  of  distributions  for  density  estimation  is 
clearly  defined,  effectively  eliminating  any  ambiguity  in  defining  the  fitting  procedure.  Further, 
unlike  the  coefficients  in  the  cubic  spline  approach,  the  various  parameters  of  the  Johnson  distribu¬ 
tion  have  specific  statistical  meaning.  Lastly,  the  Johnson  distribution  only  requires  the  estimation 
of  four  parameters  and  the  resulting  expression  can  be  easily  applied  in  equations  (4)-(8). 

In  general,  the  probability  density  function  P(r)  based  on  the  Johnson  distributions  has  the 
form 


P(r) 


7  +  6f  ( 


r~l 

A 


1  2 


(13) 


where  5  and  7  are  shape  parameters,  A  is  a  scale  parameter,  £  is  a  location  parameter,  the  function 
/(•)  is  defined  according  to  the  applied  distribution  family,  and  /'(•)  is  the  first  derivative  of  /(•) 
with  respect  to  r  [31,35].  The  distribution  families  defining  /(•),  /'(•),  and  /"(•)  are 


f(y) 


In  (y), 

In  (y  +  vV  +  !), 

y, 


for  the  Sl  (lognormal)  family, 
for  the  Sjj  (unbounded)  family, 
for  the  Sb  (bounded)  family, 
for  the  S' at  (normal)  family; 


(14) 


f'{y) 


df(y) 

dy 


< 


1  /y. 

for  the  Sl 

1  /vV2  +  !, 

for  the  Sjj 

1/[y(1  -y)], 

for  the  Sb 

1, 

for  the  Sn 

(lognormal)  family  and  y  >  0, 
(unbounded)  family  and  all  real  y, 
(bounded)  family  and  y  G  (0, 1), 
(normal)  family  and  all  real  y; 


(15) 
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< 


(16) 


f"{y) 


d2f[y) 
dy 2 


-l  A2, 

-y/(y2  + 1)3/2, 

-(l-2y)/[y(l-y)]2, 

0, 


for  the  Sl  family  and  y  >  0, 
for  the  Sjj  family  and  all  real  y. 
for  the  Sb  family  and  y  €  (0, 1), 
for  the  .S/v  family  and  all  real  y. 


The  bounds  on  these  distributions  are 


[£,  +oo) 

(— oo,  +oo) 
[£)  £  +  A] 

(— oo,  +oo) 


for  the  Sl  (lognormal)  family, 
for  the  Sjj  (unbounded)  family, 
for  the  Sb  (bounded)  family, 
for  the  Sjv  (normal)  family. 


Using  equations  (13) — (IT),  the  function  G(r)  can  be  written  as 


G'(r) 


/'(^)]2[7  +  d7(^) 


A/'  ^ 


(17) 


(18) 


When  this  expression  is  substituted  into  equation  (7),  the  resulting  stiffness  expression  as  a  function 
of  distortion  is  given  by 


[/*]  = 


r0vkT 

(/"(’-•“-«)  6 

2 

0 

s- 

+ 

C- 

CO 

S3 

1 

S3 

to 

_ ✓ 

xr  («y) 

-  a"3/2 


//  /  rQa 


-l/2_£ 
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r0a  1/2-^ 


i  2  r 


l  +  df 


v0ol 


"1/2-£ 


A/' 


t0ol 


-i/2_g 


>•  (19) 


Of  the  four  Johnson  family  distributions  initially  considered  (unbounded,  bounded,  lognormal,  and 
normal),  it  is  found  that  the  lognormal  and  normal  families  are  not  sufficiently  flexible  for  the  data 
generated  from  the  simulation  model  and  both  yield  poor  statistical  fits.  Subsequently,  only  the 
unbounded  and  bounded  families  are  considered  in  detail. 
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4.  Predictions  and  Discussion 


While  a  total  of  135  to  225  repeat  units  (Figure  2)  per  ionomer  chain  is  identified  as  the  target, 
the  chain  length  is  typically  terminated  much  earlier  due  to  chain  collision  with  a  cluster.  Actual 
simulated  chain  lengths  m  never  exceeded  45.  While  early  termination  clearly  leads  to  non-physical 
predictions  of  total  chain  length,  it  is  not  believed  that  this  early  termination  presents  a  significant 
problem  in  the  development  of  an  appropriate  P(r )  expression  because  it  is  the  lengths  between 
cluster  interactions  which  dictate  the  elastic  properties  of  the  material.  To  ensure  this,  analyses 
have  been  performed  where  all  chain  lengths  are  included  (including  those  which  terminated  early) 
as  well  as  for  the  more  physical  case  where  all  initial  and  final  lengths  are  discarded  to  eliminate  end- 
effects.  Thus,  in  the  current  study,  correction  of  free-end  effects  not  only  addresses  the  reality  that 
free  ends  do  not  contribute  to  material  stiffness,  but  also  necessarily  eliminates  early  termination 
effects.  (For  discussion  of  free  end  effects,  see  for  instance,  Treloar  [22].) 

Figures  5  and  6  demonstrate  the  effect  of  deleting  the  first  and  last  r  value  for  each  chain.  From 
these  two  figures,  it  is  evident  that  removing  the  end  effects  for  each  simulated  backbone  chain 
greatly  improves  the  behavior  of  the  data  in  the  lower  tail.  Those  r  values  computed  from  chains 
with  only  one  or  two  bonds  have  been  eliminated  and,  as  a  result,  the  second  mode  in  the  lower 
tail  of  the  density  in  Figure  5  is  not  as  prominent  in  the  density  shown  in  Figure  6.  Consequently, 
correction  of  free-end  effects  has  the  added  advantage  of  greatly  improving  the  accuracy  of  the 
Johnson  fit  since  in  general  the  Johnson  distribution  is  not  sufficiently  flexible  to  quantify  bimodal 
behavior. 

5.  Predicted  Stiffness  Values 

Applying  Equation  (9),  the  matrix  stiffness  (-E)nit  =  Em)  is  estimated  to  be  120  MPa  and  115 
MPa  for  Li+  and  Na+,  respectively,  where  these  values  are  based  on  the  experimentally  measured 
stiffness  values  (Eave)  of  75  MPa  and  80  MPa  for  the  Li+  and  Na+  cases  [8].  Figures  7-14  illustrate 
the  predicted  probability  density  functions,  with  the  location  of  rQ  indicated;  corresponding  stiffness 
values  for  the  various  cases  are  provided  in  Tables  2  and  3.  The  cubic  spline  curve  of  Figure  7 
illustrates  the  impact  of  uncorrected  end  effects  in  the  current  study  via  the  second  mode  at  very 
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(a)  (b) 

Figure  5:  Estimated  PDF  (a)  and  CDF  (b)  for  the  sodium/rectangular  case  using  a  Johnson 
unbounded  density  function. 


r  r 


(a)  (b) 

Figure  6:  Estimated  PDF  (a)  and  CDF  (b)  for  the  sodium/rectangular  case  using  a  Johnson 
unbounded  density  function  with  corrected  end  effects. 


short  chain  lengths.  Per  the  discussion  of  the  previous  section,  there  is  no  physical  reason,  other 
than  uncorrected  end-effects,  to  justify  the  existence  of  this  second  mode  as  reasonable.  Therefore, 
because  of  the  aforementioned  difficulties  associated  with  fitting  the  cubic  spline  case  even  to  reliable 
distribution  data,  no  further  cubic  spline  fits  were  employed  for  the  nonphysical  cases  (uncorrected 
end-effects). 

Because  we  expect  all  considered  cases  to  have  similar  predicted  stiffness  values  (120  MPa  and 
115  MPa  for  Li+  and  Na+  respectively),  the  bounded  Johnson  distribution  is  deemed  favorable. 
This  is  because  (a)  all  stiffness  predictions  for  the  Johnson  bounded  case  are  within  a  reasonable 
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Table  2:  Initial  stiffness  predictions  [MPa]  per  statistical  methodology,  end  effects  neglected. 


Cation/Grid  Orientation 

Cubic  Spline 

Johnson  Unbounded 

Johnson  Bounded 

Sodium,  Cubic 

3.5 

3.2 

6.4 

Sodium,  HCP 

3.4 

9.8 

Lithium,  Cubic 

-1.0 

10.4 

Lithium,  HCP 

37.8 

9.5 

ale  3:  Initial  stiffness  predictions  [MPa]  per  statistical  methodology,  corrected  end  effects 

Cation/Grid  Orientation 

Cubic  Spline 

Johnson  Unbounded 

Johnson  Bounded 

Sodium,  Cubic 

-6.6 

1.0 

9.8 

Sodium,  HCP 

26.9 

9.1 

13.8 

Lithium,  Cubic 

27.8 

4.0 

13.6 

Lithium,  HCP 

25.6 

18.1 

18.3 

range  (standard  deviation  of  3.6  MPa  from  the  mean  of  10.0  MPa),  while  the  others  have  much  more 
scatter  (standard  deviation  of  12.9  MPa  from  the  mean  of  3.7  MPa  for  the  Johnson  unbounded  case 
and  standard  deviation  of  15.9  MPa  from  the  mean  of  25.6  MPa  for  the  cubic  spline  case);  (b)  the 
previously  noted  trend  for  sodium  being  marginally  softer  than  lithium  is  most  closely  observed  in 
the  Johnson  bounded  case;  and  (c)  the  Johnson  bounded  yielded  no  negative  stiffness  predictions. 

Because  stable/repeatable  results  are  of  utmost  importance  in  the  development  of  a  new  pre¬ 
dictive  methodology,  it  is  the  stability  of  the  predictions  by  the  Johnson  bounded  case  that  is 
emphasized  rather  than  the  magnitude  of  the  predictions  themselves.  Sources  of  error  potentially 
responsible  for  the  predicted  low  magnitude  in  stiffness  include  (a)  the  semi-crystalline  nature  of 
the  polymer  matrix,  which  is  currently  neglected;  (b)  while  Figures  5  and  6  suggest  that  early 
termination  introduces  an  nonphysical  bias  toward  a  second  mode  at  very  short  chain  lengths,  and 
hence,  correction  of  end-effects  also  corrects  early  termination  effects,  it  is  conceivable  that  early 
chain  termination  interferes  with  the  prediction  of  longer  chain  lengths;  (c)  the  use  of  a  rough  weight 
average  estimate  in  the  molecular  weight  calculation  methodology  (a  more  appropriate  solution  will 
require  a  statistical  analysis  of  this  distribution  as  well);  and  (d)  the  assumption  of  spherical  cluster 
shapes,  where  oblate  or  rod-like  clusters  may  be  more  physically  accurate. 
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Figure  7:  Na+  Cubic  Configuration  with  neglected  end  effects.  PDF  curves  for  the  distribution 
of  r- values  and  corresponding  root  mean  square  value,  rQ,  for  cubic  spline,  Johnson  bounded,  and 
Johnson  unbounded  distributions. 


Figure  8:  Na+  HCP  Configuration  with  neglected  end  effects.  PDF  curves  for  the  distribution  of  r- 
values  and  corresponding  root  mean  square  value,  r0,  for  Johnson  bounded  and  Johnson  unbounded 
distributions. 


Figure  9:  Li+  Cubic  Configuration  with  neglected  end  effects.  PDF  curves  for  the  distribution  of  r- 
values  and  corresponding  root  mean  square  value,  rQ,  for  Johnson  bounded  and  Johnson  unbounded 
distributions. 
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Figure  10:  Li+  HCP  Configuration  with  neglected  end  effects.  PDF  curves  for  the  distribution  of  r- 
values  and  corresponding  root  mean  square  value,  r0,  for  Johnson  bounded  and  Johnson  unbounded 
distributions. 


Figure  11:  Na+  Cubic  Configuration  with  corrected  end  effects.  PDF  curves  for  the  distribution 
of  r- values  and  corresponding  root  mean  square  value,  rG,  for  cubic  spline,  Johnson  bounded,  and 
Johnson  unbounded  distribution. 


Figure  12:  Na+  HCP  Configuration  with  corrected  end  effects.  PDF  curves  for  the  distribution 
of  r- values  and  corresponding  root  mean  square  value,  r0,  for  cubic  spline,  Johnson  bounded,  and 
Johnson  unbounded  distributions. 
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Figure  13:  Li+  Cubic  Configuration  with  corrected  end  effects.  PDF  curves  for  the  distribution 
of  r- values  and  corresponding  root  mean  square  value,  rQ,  for  cubic  spline,  Johnson  bounded,  and 
Johnson  unbounded  distributions. 
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Figure  14:  Li+  HCP  Configuration  with  corrected  end  effects.  PDF  curves  for  the  distribution  of 
?’-values  and  corresponding  root  mean  square  value,  r0,  for  cubic  spline,  Johnson  bounded,  and 
Johnson  unbounded  distributions. 
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Figure  15:  Na+  Cubic  Configuration  with  corrected  end  effects.  PDF  curves  for  the  distribution 
of  r- values  and  corresponding  root  mean  square  value,  r0,  for  minor  variations  in  cubic  spline 
distributions  with  10  and  11  knots.  PDFs  are  nearly  indistinguishable  while  stiffness  predictions 
vary  significantly. 
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The  traditional  cubic  spline  approach  is  deemed  the  most  unreliable  of  the  considered  methods 
because  of  the  ill-defined  methodology  for  fixing  the  knots.  Any  perception  of  consistency  in  stiffness 
predictions  (based  on  the  last  3  entries  for  cubic  spline  approach  in  Table  3),  is  somewhat  misguided 
as  these  results  are  to  some  degree  ‘manufactured’  by  moving  the  knots  until  an  appropriate  response 
is  predicted.  While  this  may  be  acceptable  for  well-studied  cases,  one  of  the  primary  goals  of 
modeling  is  to  predict  those  cases  which  are  not  easily  studied  experimentally.  As  demonstrated 
in  Figure  15,  by  altering  the  placement  and  number  of  knots  we  were  able  to  alter  the  predicted 
stiffness.  In  the  illustrative  example  of  15,  the  cubic  spline  approach  not  only  yields  non-physical 
predictions  (negative  stiffness)  from  a  data  set  for  which  the  Johnson  distributions  yield  positive 
predictions,  but  also  illustrates  that  minor  variation  to  the  number  and  placement  of  knots  has  non- 
negligible  impact  on  predicted  stiffness.  The  single  greatest  challenge  with  cubic  spline  is  discerning 
the  level  of  detail  which  is  physically  appropriate.  With  the  Johnson  cases,  no  qualitative  judgment 
of  this  sort  is  required. 

It  is  proposed  that  the  reason  the  Johnson  bounded  case  offers  the  most  stable  predictions  is 
because  the  r0  position  is  consistently  nearer  the  peak  of  the  distribution.  At  this  point  in  the  P(r) 
curve,  the  expression 


r0P(r0)P"{r0)  -  rD  [P'(r0)] 2  +  P(r0)P'{r0)  <  0,  (20) 

which  is  a  necessary  condition  for  reasonable  stiffness  predictions  (based  on  Equation  (8)),  is  more 
likely  to  be  satisfied  because  P"  is  more  likely  to  be  negative  (while  we  are  generally  assured  that 
P  >  0  and  P'  <  0).  Further,  the  magnitude  of  P’  is  likely  less  sensitive  to  small  changes  in  rQ  in 
this  region. 


6.  Predicted  Stress-Strain  Responses 

Figure  16  illustrates  the  predicted  stress-strain  response  for  the  bounded  cases,  with  corrected 
end-effects  and  where  the  nominal  stress  has  been  normalized  with  respect  to  Ejnit  to  facilitate 
direct  comparison  to  experimental  data.  It  is  seen  that  up  to  several  percent  strain,  the  predicted 
stress-strain  behavior  is  consistent  with  experiment  [36].  Moreover,  while  emphasis  is  being  placed 


21 


0.25 


Figure  16:  Normalized  Stress-Strain  behavior  of  experimental  data  [29]  and  model,  both  Cubic  and 
HCP  cluster  configurations,  with  Corrected  End  Effects. 

on  the  Johnson  bounded  with  corrected  end  effects  cases,  it  can  be  shown  that  within  the  range  of 
several  percent  strain,  all  of  the  Johnson  predictions  (normalized  with  respect  to  E?;mj)  yield  similar 
stress-strain  response  trends.  The  cubic  spline  cases  introduce  considerably  more  variability  and, 
in  general,  yield  normalized  stress  predictions  that  are  considerably  higher  than  the  Johnson  cases 
at  a  given  strain. 

Further,  it  is  not  appropriate  to  compare  model  and  experiment  predictions  beyond  several 
percent  strain  because  beyond  this  point  the  strain  response  is  known  to  be  effected  by  the  hydration 
technique.  Therefore,  within  the  context  of  the  traditional  use  of  RIS  (viz.,  normalized  comparisons 
only),  the  modeling  methodology  is  deemed  appropriate  and  useful  for  study  of  solvated  Nafion, 
both  in  terms  of  treatment  of  cluster  communication  concept  for  identification  of  load  bearing  chain 
lengths  and  the  statistical  analysis. 

Traditionally,  nominal  stress  and  modulus  are  normalized  with  respect  to  ukT,  from  which 
physical  trends  may  be  inferred.  Figures  16  and  17  consider  this  approach,  where  stress-strain 
end  points  correspond  to  failure  points  (known  experimentally).  Figure  17  applies  this  approach 
to  illustrate  that  correction  of  the  end  effects  leads  to  stiffer  predicted  responses  throughout  the 
stress-strain  curve.  Similarly,  Figure  18  illustrates  that  assuming  an  HCP  cluster  morphology 
yields  stiffer  stress-strain  response  than  assuming  a  Cubic  morphology.  Therefore,  because  the 
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Figure  17:  Stress  normalized  with  respect  to  ukT  versus  strain  for  bounded  distributions,  with 
neglected  and  corrected  end  effects. 


Figure  18:  Stress  normalized  with  respect  to  vkT  versus  strain,  corrected  end  effects  with  bounded 
distributions,  for  Cubic  and  HCP  configurations. 


cluster  diameter  and  volume  fractions  are  held  constant  for  the  two  cases,  it  must  be  concluded 
that  cluster  morphology  plays  an  important  role  in  material  stiffness. 

Similarly,  the  Johnson  bounded  cases  tend  to  predict  higher  stiffness  than  the  Johnson  un¬ 
bounded  cases  (for  normalization  with  respect  to  vkT).  It  can  be  shown  that  this  trend  persists 
throughout  the  stress-strain  response  as  a  whole.  Moreover,  the  stable  nature  of  the  stiffness  pre¬ 
dictions  for  the  bounded  cases  is  mirrored  in  stable  stress-strain  trends  as  illustrated  above,  while 
the  unbounded  and  cubic  spline  cases  are  much  more  likely  to  have  erratic  stress-strain  response 
(mirroring  the  erratic  stiffness  predictions). 
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strain 

Figure  19:  Na+  Cubic  configuration  with  corrected  end  effects.  Stress-strain  predictions  normalized 
with  respect  to  Ej  and  vkT  for  minor  variations  in  cubic  spline  distributions  (refer  to  Figure  15). 

Figure  19  revisits  the  cubic  spline  cases  of  Figure  14  to  illustrate  the  corresponding  predicted 
stress-strain  response  trends.  Recalling  from  Figure  15  that  subtle  changes  to  P(r)  via  knot  selection 
lead  to  non-negligible  variation  in  the  stiffness  predictions,  it  is  observed  that  the  same  continues 
to  be  true  for  stress-strain  response  trends.  In  particular,  the  variations  in  response  illustrated 
in  Figure  19,  where  the  predictions  are  normalized  with  respect  to  ukT,  is  significant  as  this 
normalization  has  often  been  used  to  infer  physical  meanings  between  similar  distribution  data 
(i.e.,  increasing  volume  fractions  of  rigid  inclusions  [26]).  It  is  of  course  presumed  that  previous 
practitioners  of  this  modeling  methodology  have  been  aware  of  this  sensitivity  when  present  in  their 
own  systems  and  were  thus  potentially  subjected  to  a  substantial  task  in  identifying  the  proper 
cubic  spline  fit.  It  is  thus  again  emphasized  that  application  of  the  Johnson  bounded  case  represents 
both  stability  and  simplification  in  RIS  predictions  as  compared  to  the  cubic  spline  case. 


7.  Johnson  Parameter  Perturbation  Results 


Perturbation  studies  yield  insight  into  the  physical  significance  of  the  various  Johnson  parameters. 
Table  4  provides  sample  results  of  the  perturbation  studies  (for  the  case  of  Lithium,  HCP,  corrected 
end-effects,  Johnson  bounded).  In  general,  it  is  found  that  variation  of  £  (location  parameter)  has 
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the  most  significant  impact;  variation  of  5  (shape/peak  parameter)  and  A  (scale/spread  parameter) 
have  weaker  but  discernable  effects;  and  variation  of  7  (shape/skewness  parameter)  has  a  variable 
effect. 


Table  4:  Perturbation  study  of  Lithium,  HCP,  Johnson  Bounded,  corrected  end  effects. 


£ 

6 

A 

7 

rQ  (4) 

Et  (MPa) 

4.17 

2.977 

248.7 

8.689 

18.055 

16.4 

5.17 

2.977 

248.7 

8.689 

19.0275 

18.1 

6.17 

2.977 

248.7 

8.689 

20.0027 

21.1 

5.17 

2.777 

248.7 

8.689 

16.6582 

17.0 

5.17 

2.977 

248.7 

8.689 

19.0275 

18.1 

5.17 

3.177 

248.7 

8.689 

21.4726 

19.6 

5.17 

2.977 

244.7 

8.689 

18.8024 

18.3 

5.17 

2.977 

248.7 

8.689 

19.0275 

18.1 

5.17 

2.977 

252.7 

8.689 

19.2525 

18.0 

5.17 

2.977 

248.7 

7.689 

24.1299 

16.4 

5.17 

2.977 

248.7 

8.689 

19.0275 

18.1 

5.17 

2.977 

248.7 

9.689 

15.2235 

19.9 

Increasing  £  always  leads  to  an  increase  in  the  predicted  stiffness.  To  interpret  the  physical 
significance  of  this,  the  analysis  must  be  considered  within  the  context  of  Equation  (19).  From  this 
equation  it  is  seen  that  the  chain  density  v  is  held  constant  while  effectively  increasing  the  typical 
chain  length.  In  order  for  this  to  happen  both  the  material  density  p  and  molecular  weight  must 
rise,  meaning  that  any  given  bond  within  the  higher  £  material  is  carrying  less  load  for  a  given 
external  load.  This  corresponds  to  an  increase  in  the  load  carrying  capability  for  the  bulk  material, 
or  higher  stiffness. 

Increasing  5  or  decreasing  A  have  similar  effects  on  the  distribution,  where  both  lead  to  a  more 
peaked  distribution  and  marginally  higher  predicted  stiffness.  To  say  that  a  distribution  is  more 
peaked  also  suggests  that  the  molecular  weight  distribution  is  narrower.  It  has  been  found  that  a 
narrower  molecular  weight  distribution  corresponds  to  a  stiffer  material  in  other  polymers  such  as 
HDPE,  so  it  is  within  reason  that  the  same  would  be  true  for  this  class  of  polymers. 
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8.  Conclusions 


A  new  multiscale  modeling  technique  for  the  prediction  of  ionomer  stiffness  has  been  presented 
through  extension  of  traditional  RIS  theory.  Prediction  stability  is  introduced  by  fitting  the  chain 
length  distribution  P(r)  using  a  Johnson  bounded  distribution  rather  than  the  commonly  applied 
cubic  spline  approach.  For  all  considered  cases,  the  Johnson  bounded  cases  yielded  the  most  stable 
stiffness  and  stress-strain  predictions,  and  the  cubic  spline  approach  yielded  the  most  unreliable 
(variable)  results.  A  Johnson  distribution  is  further  preferred  because  its  parameters  can  be  related 
to  physical  characteristics  in  the  polymer,  whereas  cubic  spline  parameters  cannot.  Further,  the 
normalized  stress-strain  response  measured  experimentally  is  consistent  with  that  predicted  by 
the  model  up  to  several  percent  strain  (when  appropriately  normalized).  Direct  prediction  of 
stiffness  and  stress-strain  response  may  ultimately  be  possible  once  sources  of  error  such  as  early 
termination  and  estimation  of  the  molecular  weight  are  addressed  in  detail.  Lastly,  the  current 
model  gives  insight  into  significant  material  parameters.  For  example,  correction  of  end  effects 
yields  higher  predicted  stiffness;  thus  end  effects  are  non-negligible  in  this  relatively  short  chained 
system.  Further,  the  stiffness  predictions  are  sensitive  to  the  assumed  cluster  morphology,  offering 
a  tool  for  studying  this  unresolved  physical  issue. 
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